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Abstract 

We present a generalized adaptive time-dependent density matrix renormalization group 
(DMRG) scheme, called the double time window targeting (DTWT) technique, which gives ac- 
curate results with nominal computational resources, within reasonable computational time. This 
procedure originates from the amalgamation of the features of pace keeping DMRG algorithm, first 
proposed by Luo et. al, [Phys.Rev. Lett. 91, 049701 (2003)], and the time-step targeting (TST) 
algorithm by Feiguin and White [Phys. Rev. B 72, 020404 (2005)]. Using the DTWT technique, 
we study the phenomena of spin-charge separation in conjugated polymers (materials for molecular 
electronics and spintronics), which have long-range electron-electron interactions and belong to the 
class of strongly correlated low-dimensional many-body systems. The issue of real time dynamics 
within the Pariser-Parr-Pople (PPP) model which includes long-range electron correlations has not 
been addressed in the literature so far. The present study on PPP chains has revealed that, (i) 
long-range electron correlations enable both the charge and spin degree of freedom of the electron, 
to propagate faster in the PPP model compared to Hubbard model, (ii) for standard parameters 
of the PPP model as applied to conjugated polymers, the charge velocity is almost twice that of 
the spin velocity and, (iii) the simplistic interpretation of long-range correlations by merely renor- 
malizing the U value of the Hubbard model fails to explain the dynamics of doped holes / electrons 
in the PPP model. 



PACS numbers: 73. 63 Nm 



I. INTRODUCTION 



Low-dimensional strongly correlated many-body systems have been in the fore-front of 
both theoretical and experimental researches for many decades, owing to their physics being 
very different from that of the three-dimensional systems. For example, these materials 
exhibit the phenomena of spin-charge separation, wherein the spin and charge degrees of 
freedom of the electron get decoupled from each other and propagate independently with 
different velocities. From an application point of view also, these materials are in enor- 
mous demand. Amongs the genre of low-dimensional strongly correlated materials, the 
TT-conjugated polymers have attracted huge interest, being potential candidates for various 
molecular electronic and spintronic applications; examples include the organic light emit- 
ting diodes, organic semiconductors and organic thin-film transistorsii^i^ii^. However, spin 
and charge transport in these materials is still far from well understood because of the strong 
long-range electron-electron correlations that exist in these systems. Transport is strictly a 
non-equilibrium phenomena, to understand which, one needs to investigate time evolution of 
appropriate wave packets in these strongly correlated low-dimensional systems. The density 
matrix renormalization group (DMRG) technique advanced by Steve White^i^, has proved to 
be a very powerful numerical method for studying large interacting low-dimensional quantum 
lattice systems. Originally formulated as a ground state method, this technique has been 
mostly used for studying static (equilibrium) quantum many-particle phenomena. Later, it 
was extended to calculate frequency- dependent spectral functions by the correction vector or 
the Lanczos techniques^'^'^'^'^'^'^,— , thereby adapting it to deal with dynamical (equilib- 
rium) many-body phenomena. In this regard, the dynamical DMRG (DDMRG) method^ 
turned out to be the spectral method, best suited for obtaining extremely accurate spectra. 
However, the DDMRG technique is limited to only one momentum and a narrow frequency 
range at a time. Constructing an entire spectrum using this method in order to obtain a 
reasonable grid in frequency and momentum space involves independent runs for each fre- 
quency and momentum and is therefore, computationally highly intensive. An alternative 
route exists for obtaining an entire spectrum in a single run which involves time evolving an 
appropriate wave packet in real space and time, followed by converting the information from 
(r,t) space to {k,uj) space using double Fourier transform. Quantum dynamics of a wave 



packet in real space time is needed to obtain the {f,t) data, which for large systems was not 
feasible, until recently. Three of the recent time-dependent DMRG (t-DMRG) techniques 
are the pace keeping DMRG scheme due to Luo, Xiang and Wang (LXW)^, the adap- 
tive DMRG (t-DMRG) method^»^>^ and, the time-step targeting (TST) technique due to 
Whitens. In this paper we present a time- dependent DMRG scheme which is a hybrid of the 
LXW and TST algorithms. The organization of the paper is as follows: In Sec. I we give an 
overview of the existing time-dependent DMRG techniques, followed by a detailed discus- 
sion of the LXW and TST algorithms in order to compare their strengths and weaknesses. 
Section II is devoted to a detailed presentation of our algorithm. Section III discusses some 
numerical issues associated with the DTWT scheme. In Sec. IV we compare the technique 
with the LXW and TST schemes. In Sec. V we present real time dynamics of spin-charge 
separation in regular Pariser-Parr-Pople (PPP) chains using our DTWT algorithm. Section 
VI provides the summary and conclusions of this work. 

A. Overview of existing time-dependent DMRG methods 

The study of out-of-equilibrium phenomena in strongly correlated low-dimensional sys- 
tems has attained a forefront; DMRG has played a significant role in this context. Real 
space-time quantum dynamics using DMRG was introduced by Cazalilla and Marston^. 
They calculated time evolution of one-dimensional systems under the influence of an ap- 
plied bias. Their approach involved the use of DMRG for constructing the Hilbert space of 
the Hamiltonian using only the ground state, and subsequently solving the time-dependent 
Schrodinger equation numerically using the Hamiltonian matrix obtained in a fixed basis. 
Hence, the approach of Cazalilla and Marston is essentially static with respect to the Hilbert 
space in which time evolution is carried out. It is expected that when the evolving wave 
function becomes significantly different from the ground state, i.e., "it moves out of the 
Hilbert space used for time evolution," it will loose accuracy. However, in the systems stud- 
ied by them, this (expected) loss in accuracy with time did not occur within the time period 
for which they carried out the time evolution. The drawback of performing time evolution 
using a fixed DMRG basis is, one needs to keep a substantially large number of the density 
matrix eigenvectors (DMEVs), m, so that the evolved state is well described by the DMRG 



basis at large times. Luo, Xiang and Wang^^ showed an effective way, also known as pace 
keeping DMRG, to construct DMEV basis for time evolving a wave function over a longer 
time interval. This is done by constructing a weighted average density matrix from the time 
evolved wave functions at discrete time steps Ar in the time interval 0- T. Thus, 

Nt 

PL/R = TrR/L^a, I ip{ti)){'ilj{U) I; ^a, = l, (1) 

i=Q i 

where Nt = The weights are taken to be 1/2 for i = and otherwise. However, 
their approach suffers from two significant drawbacks. The LXW scheme performs full time 
evolution of the wave function at each infinite DMRG step thereby making it extremely time 
consuming. In DMRG calculations with multiple target states, truncation error is reduced 
by keeping the number of DMEVs (m), greater than the number of target states. But in 
the LXW algorithm, the number of target states Nt [each target state corresponding to 
'4){ti) at time tj] is usually » m, thereby decreasing accuracy. Modification of the original 
LXW scheme in which Pl/r is obtained as an average in each sub- interval At, where At = 
pAr, p << Nt leads to a decrease in number of target states required, and consequently 
increases accuracy of the results^. In an earlier study, we demonstrated this for T = 33 
fs, Nt = 50,000, Ar = 0.00066 fs and p = 200-500. In fact, the accuracy is not degraded 
even for p = 100. However, our earlier work has shown that this modification works well 
for nearest-neighbor Hamiltonians like the tight-binding (Hiickel) and Hubbard models only. 
But, in the case of models with long-range interactions, like the PPP Hamiltonian, the above 
modification fails and we need to average the density matrix over all time steps. 

In the context of simulating time evolution of matrix product states (MPS), Vidal de- 
veloped a novel numerical scheme called the time-evolving block decimation technique^^. 
As DMRG is closely related to MPS, this method was immediately utilized by the DMRG 
community to develop a very powerful numerical technique called adaptive t-DMRG^>^, for 
studying real-time quantum dynamics in strongly correlated many-particle systems. The key 
idea of t-DMRG technique is to incorporate the Suzuki- Trotter (ST) decomposition of the 
time evolution operator exp(— zAriJ) into the finite DMRG algorithm. Usually, second-order 
ST decomposition is used but higher order ST decompositions of the matrix exponential can 
be employed as well. The second-order Suzuki- Trotter breakup is 



where Ha is the Hamihonian for the bond connecting sites (2j-l, 2j) and Hb connects sites 
(2j, 2j+l); j = 1,2 ■■■ for Ha and j =1,2, ■■■ y-l, for Hb- This breakup leads to 
grouping of terms such that all terms within Ha or within Hb which commute with each 
other, leading to high accuracy in the Trotter breakup. A DMRG super-block state, at a 
particular step n of the finite-system sweep, can be represented as 

L R 

where | L) and | R) represent the truncated DMEV basis states of the left and right blocks, 
and I represent the Fock-space states for the two central sites. Thus, an operator 

(A) acting on the two central sites can be exactly expressed in the same optimal basis as 

The ST decomposed time evolution operator for the bond connecting sites n and n+1 can 
be applied exactly on the DMRG super-block state at step n of the finite DMRG algorithm. 
A full finite DMRG sweep corresponds to a time evolution by 2Ar of the full system. In 
the t-DMRG scheme, two types of errors are involved: (1) the DMRG truncation error 
and (2) the ST decomposition error, which for second-order decomposition is {OAt^), in 
each time step. Adaptive t-DMRG based on the Suzuki-Trotter decompositions of the 
time evolution operator is restricted to nearest-neighbor interactions only since the break- 
up into Ha and Hb is valid only in this case. Schmitteckert^ proposed a Krylov space 
approach for obtaining the exponential time evolution operator. In this method a small 
matrix representation of the original Hamiltonian, H, is obtained in the basis vectors of 
the Krylov space, namely {"a/", Hl^, H^lt, ■ ■ ■ , H''~^lt] "a/", an arbitrary vector} for an Ith 
dimensional Krylov space. This permits explicit form of the time evolution operator in a 
finite optimized basis. He, however, employed the LXW scheme for obtaining the DMRG 
space spanned by the Hamiltonian. 

White proposed a second adaptive approach called the time-step targeting (TST) tech- 
nique to circumvent the nearest-neighbor limitation associated with ST break-up^. Two 
key ideas are introduced in this scheme. The first being, there is no conceptual need to 
have the same time step for time evolution and for building the Hilbert space of the evolv- 
ing wave function. This simply implies that the time-scales Ar on which time evolution 



is discretized and At on which the DMEV basis is adapted, need not be identical (this is 
the crucial difference between the LXW and TST schemes). The observation from the work 
of Cazalilla and Marston that even a static state space remains a good choice for some 
finite time, implies the choice of At > Ar. Thus, one can evolve the wave function using 
any of the non- Trotter methods, the one used by White and Feiguin being fourth-order 
Runge-Kutta (R-K) technique. Other available methods are the Crank-Nicholson scheme or 
Krylov-space-based decomposition of the matrix exponential^. The second key idea of the 
TST technique is to adaptively build the Hilbert space representing the time evolving wave 
function, at intervals At, from the states which are expected to appear in future through 
time evolution. For this purpose, several DMRG sweeps are carried out at a fixed time t. 
For each (L — — an+i — R) configuration during these finite sweeps, a fourth-order Runge- 
Kutta integration is carried out, and the time-evolved states | ip(t + jAt)), j = 0, |, |, 1 
are used to build the reduced density matrices of the blocks using Eq. (1). After sufficient 
number of sweeps, when the Hilbert space optimally represents the wave function | ipit)), 
final time evolution is performed using a time step Ar ^ The new wave function is 
then used to build the DMEV basis for the next time propagation. Both LXW and the 
TST techniques are generalized adaptive time-dependent DMRG schemes since they can be 
applied to any Hamiltonian, with arbitrary range (beyond nearest neighbor) of interactions. 
This flexibility stems from the fact that the time evolution operator is not decomposed using 
Suzuki- Trotter-type decomposition schemes. In the LXW technique, adaptive construction 
of the Hilbert space of the time-evolving initial state as well as temporal propagation are 
performed within the context of infinite DMRG algorithm. On the other hand, TST tech- 
nique uses the finite DMRG scheme to adaptively update the Hilbert space of the time 
evolving wave function as well as, propagate it in time. The time steps used for state space 
updating and time evolution of the wave packet are the same in LXW scheme while they 
are chosen to be independent in the TST scheme. As stated by White, the time step for 
basis adaptation (At) is larger than that for evolution (Ar). 

We have compared the two techniques and have found that, for a given value of m and 
system size N, LXW scheme besides being comparable in speed with the TST algorithm, is 
also more accurate. In the LXW scheme, at every system size we need to evolve the wave 
packet over the entire time period, before getting the DMEV basis for moving to the next 



system size. In the TST scheme, although we do the time evolution only for the desired 
finite-system size, we employ a finite DMRG scheme at every time evolution step At, which 
entails large CPU times. The higher accuracy of the LXW scheme arises from the fact 
that DMEV basis is obtained from a weighted average density matrix, constructed from 
time-evolved wave packets at all time intervals, while in TST scheme the weighted average 
density matrix is constructed from the wave packets 0, and At for each time step 

of evolution At. 

If one uses the window modification of the LXW technique as proposed by us2^, thereby 
decreasing the number of target states, the accuracy as well as computational speed can be 
further improved. The above observations leads one to conclude that the LXW technique is 
superior compared to the TST scheme. However, LXW scheme also suffers from two serious 
drawbacks: first, LXW method is not quasiexact unlike the TST scheme. An approximate 
DMRG algorithm is quasiexact when the error in the observables is strictly controlled by the 
truncation error, e = 1.0 — XliP*' Pi being the dominant eigenvalues corresponding to the 
reduced density matrix eigenvectors which are retained. In case of a quasiexact scheme, the 
errors in expectation values of the system's properties are either proportional to e or i/e, that 
is, {{0{t)) DMRG — {0{t)) exact) oc 6 or y/l. The infinite-system DMRG technique applied 
to a finite system is not a quasiexact approximation scheme even though e goes to zero as 
m increases. However, in the absence of any metastable ground states and with "sufficient" 
number of sweeps, the finite system DMRG algorithm is quasiexact. LXW scheme being an 
infinite system DMRG algorithm applied to finite systems, is nonquasiexact^S. Second, in 
case of long-range Coulomb interaction as in the PPP Hamiltonian, our window modification 
fails, thereby making it necessary to retain all the time-dependent target states. This makes 
the LXW scheme computationally inefficient. Finally, both the LXW and TST techniques 
are computationally time consuming. These observations motivated us to develop a real- 
time evolution method which will have the strengths of both these methods while overcoming 
their limitations of long computational times and poor accuracy, especially for long-range 
interacting models. 



II. DOUBLE TIME WINDOW TARGETING (DTWT) ALGORITHM 



Conceptually, the scheme which we have implemented is a hybrid of the LXW and TST 
algorithms described in the previous section. These schemes differ from each other mainly in 
the prescription for constructing the weighted average density matrix. The different schemes 
are compared schematically in Fig. 1. In the LXW scheme, the weighted average density 
matrix is obtained from the density matrices of the time evolved states at all times ( 0, T) in 
steps of At, the total time interval from ^ T being At. In the TST scheme, in each time 
interval At, a weighted average density matrix is constructed, typically with the states at the 
beginning and end of the time interval and at two equidistant intermediate points. The time 
evolution is carried out with a time step At such that. At is ~ within the interval. For 
time evolution over the next time step At to 2Ar, the density matrix is constructed using 
the time interval from At to {At + At). Besides, in the TST algorithm, unlike in the LXW 
algorithm, finite DMRG procedure is carried out at each time at which the density matrix 
computation is carried out. In the present scheme which we call the DTWT technique, we 
construct the density matrix over a time interval 2pAr, as a weighted average of the density 
matrices built at time steps of length At, and use the DMEVs from the resulting density 
matrix to constructing the Hilbert space of the Hamiltonian, for time evolving the desired 
wave packet from At to pAr. In this technique, each time interval of length At is broken 
into p sub-intervals of length At such that, while the wave packet evolves by At, the basis 
is adapted over an interval 2At. After every time evolution by At, the interval is slided 
by 2At for constructing the Hilbert space for the next time evolution. Thus, approximate 
future states for a time period 2At are used for evolving the wave packet over a time interval 
At, in steps of At. As in the TST scheme, we employ finite procedure to get accurate wave 
functions in each time interval. 

A. Initial wave-packet construction 

A conventional infinite DMRG algorithm is carried out to build the system of desired 
final size on which time evolution is wished to be performed. We wish to study the dynam- 
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FIG. 1: Pictorial representation of the construction of the time averaged reduced density matrices 
for left/right block (pl/r) in the (a) original LXW, (b) modified LXW, (c) TST, and (d) DTWT 
schemes. In all the schemes, At = and At = pAr. In the original LXW scheme (a), At 
= At {p = 1) is used while in our modified LXW scheme (b), At = pAt, with p ~ 100 is used. 
In the TST scheme (c), a sliding time window of length At is used for updating time step At 
(At ~ IOAt). In the DTWT method (d), a sliding time window of length 2 At = 2pAT, is used 
for updating time window of length At = pAt, p ~ 100. In the last two schemes, finite DMRG 
procedure is carried out to obtain basis adaptation and time evolution. 



ics of spin and charge transport in a one-dimensional model with long-range interactions. 
Therefore, the initial wave packet at t = is formed by annihilating an upspin ('|") electron 
from site 1 of a neutral system in the ground state (| ipgs)), 

I ^(0)) = ai,t I ^;,). (5) 

Since we are interested in the time evolution of this charged wave packet formed by anni- 
hilating an electron in the chain, the Hilbert space formed from the DMEVs of the neutral 
system alone, would be inappropriate. Hence we construct the half-block density matrices 
at each system size as a weighted average of the density matrices of the desired wave packet, 
the neutral ground state and other relevant states. These density matrices correspond to 
system of the same size but with different particle numbers. Thus, at each step of the infinite 
scheme, the half-block reduced density matrices are formed in the following way: 



Pl/r = Tr^/L^^wo I ^{0)){ij{0) I + I tpj){^j 1^ 



= Ttr/lIuo I m){m I + I I + E^^- I ^^■)(^^- I ' 

V j=2 ^ 

where, | ip{0)) is the initial state with weight uq, \ tpj) is the jth relevant state having 
weight Uj, and + J2]=i — ^ ^^'^ ^^e number of relevant states including the 
neutral ground state is r. Other states which can be considered as relevant target states 
are the ground state of the ion, especially in case of inhomogeneous systems. From our 
earlier studies on the LXW algorithm we have found that for systems such as (CN)^; and 
(PN)^ consisting of more than one type of atom, it is necessary to target the corresponding 
cationic (| 'ip^s)) or anionic (| ip^s)) ground states also, in addition to the initial wave packet 
and neutral ground state. Thus, in this case the density matrix is given by 

PL/R = TrR,L[uj, I m){m I + ^1 I I + <s I i^ts)(^ts)y (7) 

An exhaustive analysis of the dependence of the charge densities and spin densities of the 
initial wave packet on the weights of the target states showed that the nature of the initial 
state does is not very sensitive to the weights. Hence, as in the LXW scheme, we have 
chosen ujq = 0.8 and = u'^^ = 0.1 for all the systems we have studied. 



The initial state for the final lattice, which is obtained using the infinite DMRG scheme, 
is numerically evolved in time using the time dependent Schrodinger equation (TDSE). 



The purpose of this time evolution is to obtain a set of time- dependent states, with which 
the initial Hilbert space of the time-evolving wave packet can be constructed for starting 
the finite DMRG scheme. The TDSE has been solved numerically using the second-order 
multistep differencing scheme (MSD2)^. We settled on this scheme after trying more ac- 
curate schemes such as fourth- and sixth-order multistep differencing schemes (MSD4 and 
MSD6),— , and the fourth-order R-K method. The final evolution step is the one that needs 
accurate integration of the TDSE. The MSD2 scheme is given by the following equation: 

I ij{t + mAT)) = -2iH I tlj{t+[m-l]Ar))+ \ tlj{t+[m-2]AT)) +0{AT^);m e [2,2p] (8) 

The Hamiltonians used in the present study are time-independent and H = {H~^ — Eq), 
where is the Hamiltonian of the positively charged system (cation) and Eq is the ground 
state eigenvalue of H^. The flow-chart for the initial infinite DMRG procedure is given in 
Fig. 2 

B. Hilbert-space adaptation and time evolution using finite DMRG scheme 

After obtaining the initial wave packet, the neutral and the ionic ground states, and the 
time evolved wave packets | ip{kAT)), k G [l,2p], using infinite DMRG algorithm (see Fig. 
2), for a system of desired size, finite DMRG procedure is carried out to adaptively construct 
the Hilbert space of the Hamiltonian and propagate the wave function in time. Each finite 
DMRG step generates a system configuration having two blocks joined together by two single 
sites. Two full sweeps of the finite DMRG procedure involves 4(^ — 2) steps, Ns being the 
final system size; we call (^ - 2) finite DMRG steps a half sweep. Hence, two full sweeps 
imply sweeping from (^ - 1, 1, 1, ^ - 1) ^ ■ • ■ (iV, - 3, 1, 1, 1) ^ ■ ■ ■ (^ - 1, 1,1,^- 1) 
{1,1,1, Ns — 3) — 7- ■■■ (^ — 1,1,1,^ — 1), where the first and last numbers in 
parenthesis give the number of sites in the left and right blocks, respectively. At every step 
of the finite DMRG procedure starting from infinite DMRG solution corresponding to the 
system (^ — 1,1,1,^ — 1), the following operations are performed: 

(1) At time t = 0, the reduced density matrix for the appropriate (left or right) block 
is computed using the neutral ground state | ?/;° ) with an assigned weight , the initial 
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FIG. 2: Flow-chart for the infinite DMRG algorithm in the DTWT technique, k = 1,2,3, , 

(^-1), is the half-block length; Ng is the total number of sites in the final lattice, a and a' 
are the two new sites added during the infinite DMRG scheme; ipgsii^ + + 1)'] ^i^d 

ip'^s[{k + l),a,a' , {k + 1)'] are the neutral and ionic ground states, ^[(/c + l),a,a' , {k + l)';mAr] 



is the time evolved wave packet at i = mAr, obtained using the MSD2 technique. 



wave packet | ip{0)) with weight uq, the ionic ground state | ip^^) with weight and all 
other preliminary time evolved states | ip{kAr)) [k = 1,2,- ■ ■ ,2p) with weights Uk- The total 
weight is normalized to unity [Eq. 



PL/R = Trn/L[uj% I 0«) + | ^+)(^+ | + ujo \ m){m I 

X (9) 
+ ^Wfc I + A;Ar))(^/'(t + A;Ar) I j. 
k=i ' 

At other times it 7^ 0), the reduced density matrix for the appropriate (left or right) block is 

constructed using the initial state, the ionic ground state, and 2pAr time-dependent states 

(the neutral ground state is not considered), as given below: 

PLIR = Trn/L^uj^s I i^t)('p9s I + c^o I ^(0))(^(0) | 

(10) 



+ J2^k\'ip{t + kAT)){'4}{t + kAT)\\. 
k=l ' 



The weights w^^, w^j,, and are adjusted (as in the infinite DMRG scheme used to generate | 

tlj{t + kAt)) , k = 1,2,3, ,2p, described in the previous section), after exhaustive analysis. 

In the case of TST algorithm, the weights of the target states can all be chosen to be equal or 
unequal. However, in our case we find that equal weightage of all the target states severely 
deteriorates the results. Our tests have shown that the optimal unequal weights for Eq. (9) 
are wn = 0.7, and = ojj; = 0.1 so that Uk = — — — — — — = — , V k. In case of 

^ ^ gs gs t p p ' 

Eq. (10), the optimal weights are coq = 0.8, uj+ = 0.1, and Uk = ^^'^ ' ~ "^''^ = ^, V A;. 

(2) The resulting left- or right-block density matrix is diagonalized and the DMEVs 
corresponding to the double time window, t ^ (t + 2pAr) are chosen. Using this DMEVs, 
the renormalized block Hamiltonians (ff^/^) and operators [Ol/ji) are constructed. Using 
these renormalized Hamiltonians and operators, the super-block Hamiltonian for the current 
system configuration with required particle numbers is obtained. 

(3) The super-block Hamiltonians for neutral and ionic systems are diagonalized us- 
ing either the Lanczos^, Davidson's^^, or any other iterative sparse matrix diagonalization 
algorithn>2i to obtain the ionic and neutral ground states. In case of the time evolution 
from t = to t = At, both the neutral and the ionic ground states are needed while for the 
subsequent time evolutions {t 7^ 0), only the ionic ground state is required. This implies 



that two super-block diagonalizations are needed for the initial time evolution while only 
one super-block diagonalization for the remaining evolutions, at every finite DMRG step. 
We have used the Davidson's algorithm for diagonalization of the super-block Hamiltonian. 

(4) For time evolution from t = to t = At = pAr, the initial wave packet is explicitly 
constructed from the neutral ground state according to Eq. (5). Hence, the neutral ground 
state is retained as a target state in Eq. (9). For the remaining time evolutions {t 7^ 0), 
the initial state is transformed from the old DMEV basis to the current one using White's 
wave function transformation^; therefore the neutral ground state is no longer required. 

(5) The initial wave packet and super-block Hamiltonian which are both expressed in the 
current DMEV basis, are used to solve the TDSE numerically using the MSD2 scheme from 
t — (t + 2At) in time steps of At. 



(6) After two full sweeps the final configuration (^ — 1,1,1,^ — 1) is evolved from t — )■ 
(t + At) = {t + pAr). This is the final (single) time window evolution. This is carried out 
using the fourth-order R-K method. The fourth-order R-K method propagates | ip{t)) to 
I ip(t + At)) using four vectors | /ci), | ^2), | k^), and | /C4) defined as 

I ki) = -lATHit) I tl){t)), 
\k^) = -iATH{t + ^)(^\ij{t)) + l/2|A;i)^, 
\k,) = -^ATHit + ^)(^\m) + 1/2 I A;^)), 

\k^) = -iATH{t + AT)(^\ij{t)) + |A;3)^. (11) 

The state at time (t+Ar) is given by 

I ^{t + At)) ^ i (^1 fci) + 2 I A;2) + 2 I ^3)+ | k,)^ + 0{At'). (12) 

The obtained time-dependent states can be either saved to disk for future use or they can 
be used on the fly for calculating dynamical observables. The final state | ip{t + pAr)) is 
then used as the next initial state and the same procedure is repeated for the next single 
window time propagation (see Fig. 3). 



In the infinite DMRG part of the DTWT technique, the half-block density matrices and 



the corresponding DMEVs are constructed from the initial wave packet, the neutral ground 
state and the ionic ground state, all of which are time-independent real wave functions. 
Hence in the infinite DMRG part of our algorithm, all operations are performed in real 
arithmetic and the quantities (scalars, vectors, and matrices) involved, are all real. However, 
in the finite DMRG part we encounter time-dependent states which are complex quantities, 
resulting in complex reduced density matrices, DMEV basis, renormalized left- and right- 
block Hamiltonians, block operators, and super-block Hamiltonian matrices. Therefore, the 
DTWT method as implemented by us performs real arithmetic for the infinite DMRG part 
and complex arithmetic for the finite DMRG part. 
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FIG. 3: Basic scheme showing the use of finite DMRG algorithm in the DTWT technique. Q + k 
{Q' — k) is the length of left (right) block. rjF denotes the direction of sweep: rjp = 1 (-1) implies 
left (right) — )■ right (left) sweep, hlfswp = number of half sweeps. Time evolution is performed 
using fourth-order R-K technique while basis adaptation is done using the MSD2 technique. 



III. SOME NUMERIAL ISSUES 



It is well known in the literature that the sparse super-block matrix diagonalization is 
the most time-consuming step in both the infinite and finite DMRG schemes. For the time 
evolution from t = to t = pAr, each step of the finite part of our DTWT algorithm 
involves two sparse matrix diagonalizations: one for the neutral system's Hamiltonian and 
the other for the ionic system's Hamiltonian. For all subsequent time evolutions, t > 
pAr, diagonalization of only the super-block Hamiltonian for the ionic system is required 
at each step of the finite part of the DTWT procedure. If the sparse matrix diagonaliza- 
tions involved at each step of the finite DMRG part of our algorithm can be replaced by 
the wavefunction transformation introduced by White^, computational time is substantially 
reduced. However, our tests indicated that for the finite DMRG steps, [Ng — 4, 1, 1,2) ^ 
{Ns — 3, 1, 1, 1) and (2, 1, 1, A'^s — 4) ^ (1, 1, l,Ns — 3), the wave function transformation fails 
to be accurate. Thus, the matrix diagonalization step needs to be selectively replaced by 
the wave function transformation. Taking recourse to White's wavefunction transformation, 
we reduce the number of sparse matrix diagonalizations for the initial time evolution {0 — t- 
pAr), from 8(^-2) to 8 for two full sweeps; for all subsequent time evolutions (t > pAr), 
the wave function transformation reduces the number of super-block Hamiltonian diagonal- 
izations from 4(^-2) to 4, for (every) two full sweeps, thereby reducing the computational 
time significantly. Obtaining the ground state wave function and energy through matrix 
diagonalization involves the solution of an eigenvalue equation, \ ipgs) = Egg \ ipgs)^ 
while obtaining the same using the wave function transformation involves transformation of 
the ground state from the old to current DMEV basis, followed by the evaluation of Egg as, 
Egs = ■^^^r^/T^i ^ being the super-block Hamiltonian in the current DMEV basis. 

As the number of target states used for constructing the reduced block density matrices 
increases, the overall computational time for the infinite and finite DMRG schemes also 
increases. The infinite part of DTWT algorithm uses three target states while the finite part 
needs 2p states as target states. Hence, the finite part of our algorithm is slower compared 
to the infinite part. To overcome this problem, we have used the window modification that 
we developed and employed in the context of LXW algorithm. Instead of retaining all the 
2p time-dependent states as target states, we keep ^ target states; instead of 10, other 



variations can also be employed. Incorporation of the window modification into the finite 
part of our algorithm has an interesting consequence, namely, the 2p < m condition can be 
replaced by < m, thereby making it feasible to increase the length of both the double 
(2At) and single (At) time windows. 

The time step error associated with the fourth-order R-K technique is O(At^) while that 
associated with the MSD2 scheme is 0{At^). Thus, the former time propagation scheme 
is more accurate than the latter. However, the TDSE solver required for constructing the 
Hilbert space of the time-evolving wave function need not be very accurate compared to 
TDSE solver needed to propagate the wave packet in time. Hence we have used two different 
time propagation techniques for the two different windows present in our algorithm. For the 
double time window (2At) evolution which is used for basis adaptation, we have used the 
second-order multistep differencing (MSD2) scheme^^, given by 

\ilj{t + AT)) = -2tHAT \ip{t))+ \ip{t- At)), (13) 

which uses one sparse matrix-vector multiplication (SMVM) operation for every propagated 
time step At. The two states | ip{t)) and | tlj{t — At)) are obtained once in the beginning of 
every time evolution using fourth-order R-K technique. Hence, every time evolution from t 
— 7- (t+2pAr) involves (8 + [2p — 2]) = {2p — 6) SMVMs. For the final single time window 
propagation we employ the fourth-order R-K technique as the TDSE solver. Hence, every 
time evolution from t — (t + pAt) would require 4p SMVMs. Therefore, each single time 
window propagation using the MSD2 and R-K techniques involves ^[2j9 — 6] x — 2] +4p^ 
SMVMs. 

IV. COMPARISON OF THE DTWT SCHEME WITH THE LXW AND TST 
TECHNIQUES 

In this section we present a comparative analysis of the computational efficiency of our 
DTWT scheme with the LXW and TST schemes. In the TST technique, as pointed by 
White, the wave function is propagated in time after either one or several half-sweeps de- 
pending on the how many half sweeps are needed to update the Hilbert space of the wave 



function adequately. Each half sweep (according to our definition of half sweep) needs 
— 2) finite DMRG steps. If the wave function is propagated in time by At after every 
half sweep, then for the total time propagation, Np(^ — 2) finite DMRG steps are needed. 
However, if the wave function is propagated after every q half sweeps, then Npq{^ — 2) 
finite DMRG steps are needed for propagating the initial wave packet over the entire time 
interval, to T. In the DTWT technique, for every two full sweeps {q = 4) the wave function 
is evolved by a time interval At which implies that Nq{^ — 2) DMRG steps are needed for 
evolving the initial state over the entire time interval. Therefore, the ratio of DMRG steps 
in DTWT to TST is either ^, if g 7^ 1, or (4/p), if g = 1 in the TST procedure. Since 
p is chosen to be > 10^, the DTWT scheme is more efficient than the TST technique by a 
factor of > 25. The actual CPU time involved in the time evolution by the fourth-order R-K 
procedure is the same in both TST and DTWT schemes. In case of the LXW algorithm, 
the DMRG steps involved are fewer for a system size Ng, the number of DMRG steps is 
only — 1). However, we evolve the wave packet at each system size and this is the CPU 
intensive part of the calculation. Besides, the target states in LXW scheme is equal to num- 
ber of time intervals (specially for systems with long-range interactions), which is ~ 10^ 
or more. Therefore the DMEV cut-off required for comparable accuracy is huge (~ 10"*) 
and unattainable. This leads to large errors in the LXW scheme at long times. To test the 
method we compared the DTWT results for various m values with exact time evolution of a 
14-site Hiickel as well as PPP chain. We also compared these results with those obtained by 
the TST method for a 14-site Hiickel chain and LXW for a 40-site Hubbard chain with = 
4. To test whether our DTWT scheme is applicable to other geometries, we have compared 
the exact time evolution of a 12-site biphenyl molecule and a 14-site stilbene molecule (see 
Fig. 4), both modeled by the PPP Hamiltonian, with different DMEV cut-offs, m. We have 
dealt with Hiickel chains of length up to 14-sites since the number of states (although known 
exactly) become too large and the computations become cumbersome for time propagation 
of the wave packet. 

The model Hamiltonians used in this study are the tight-binding Hamiltonian^i^, also 
known as the Hiickel Hamiltonian to chemists, the Hubbard^^i^i^^ and PPP Hamiltonian^^i^. 



The second quantized Hamiltonians for these models are give below^: 

-f^Hiickel = ^ ] ^ij {^l,(T^j,o^ ~^ ^],a^i,'^)^ 

<«i>,o-=t,^ 

-f^Hubbard = -f^Hiickcl + UiUi^^rii^i, 

i 

iJppp = i^Hubbard + ^Vij{ni- Zi){nj - Zj). (14) 

i>j 

In the first equation, the summation is restricted to bonded neighbors < ij >, tij is the hop- 
ping integral between bonded neighbors. Ui is the on-site Coulomb repulsion term (Hubbard 
f/-term) of the ith site, c]^ (Q,cr) creates (annihilates) an electron with spin a at (from) the 
ith site, the corresponding number operators, rij = (rij^o- + rii^^cr) is the charge 

density at site i, Vij is the intersite Coulomb repulsion between lattice sites and Zi 

is the on-site chemical potential of site i. For homogeneous systems, f/j = U, and is a 
measure of the Coulomb repulsion between two electrons of opposite spins occupying the 
same site. [U /t) characterizes the strength of electron correlations. In case of homogeneous 
sp^ carbon systems, all sites are singly occupied for charge neutrality and hence Zi = l,\f i. 
The inter-site electron-electron repulsion term (Vij) in the PPP model is phenomenologically 

2 

interpolated between U for zero separation and f- for inter-site separation rij — )■ oo; thus, 
the explicit evaluation of the repulsion integrals is avoided. There are two widely used inter- 
polation schemes for evaluating Vy-. Ohno scheme^ and, the Mataga-Nishimoto scheme^. 
In the Ohno interpolation scheme which we use, the Vij term is given by 



Vij = 14.397 



28.794 



2 -| -1/2 

+ r^- 



Ui + Uj, 

and decays more rapidly than the Mataga-Nishimoto formula which is shown below. 

2.0 r,; 



(15) 



Ui + Uj ^ 14.397 



(16) 



In both the above interpolation formulas, rij is measured in angstrom (A) while U and Vij 
are measured in electron volt (eV). In the Hubbard model calculations, we have used = 

4 while in the PPP model we have used standard parameters {U = 11.26 eV, t = -2.4 eV, 
and r = 1.397 A; r being the uniform C-C bond length) and Ohno parametrization. In Fig. 

5 we compare the DTWT and TST schemes for different DMEV cut-offs m for a Hiickel 
chain of 14 sites, with exact results. From Fig. 5 it is clear that for a given system size and 
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FIG. 4: Molecular structures of (i) biphenyl and (ii) stilbene molecules. Biphenyl and stilbene 
are the monomers of polyparaphenylene and polyparaphenylvinylene polymers, respectively. Site 
numbering scheme as employed by us, is also shown in the figure. 

m, the DTWT algorithm has a better accuracy than the TST algorithm. The results for 
the Hubbard chain (Fig. 6) of 40 sites show a smooth convergence in the DTWT scheme 
as m is increased. This gives confidence in the DTWT scheme. The LXW method with 
m = 200 differs from our results quantitatively at long times (Fig. 6). Our results for the 
14-site PPP chain indicate that long-range interacting models require higher m for attaining 
the same convergence as with the nearest-neighbor models (Fig. 7). In Fig. 8 we show the 
time evolution of the charge and spin densities at the sites numbered according to Fig. 4. 
Fig. 8 also shows that molecular topologies which are not linear, require larger DMRG basis 
dimension to attain accuracies comparable to open chain system. 
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FIG. 5: (Color online) Time evolution profiles of charge (left curve) and spin (right curve) densities 
at first and last sites of a 14-site Hiickel chain, using three different techniques. The color codings 
for the curves are as follows: black curve with circles = exact; red curve with squares=TST(m=64); 
green curve with left triangles=TST(m=100); blue curve with right triangles=DTWT(m=64); or- 
ange curve with up triangles=DTWT (m = 100). The m=100 DTWT curve is almost indistin- 
guishable from the exact curve. At = 0.066fs and At = ^ are used in the TST technique. In the 
DTWT scheme, At = 0.066fs is employed. 
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FIG. 6: Comparison of time evolution of charge (left curve) and spin (right) densities at the first 
and last sites of a 40-site Hubbard chain, using the LXW and DTWT techniques. In both the 
DTWT and LXW schemes. At is taken to be 0.066fs. 
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FIG. 7: (Color online) Comparison of exact versus DTWT time evolution of charge (top curve) 
and spin (bottom curve) densities at first and last sites of a 14-site PPP chain. Ar is chosen 
to be 0.0066fs. The color coding is as follows: black curve with circles = exact; red curve with 
squares=DTWT(m=64); green curve with up triangles=DTWT(m=100); blue curve with down 
triangles=DTWT(m=150); orange curve with left triangles=DTWT(m=200); violet curve with 
right-triangles=DTWT(m=250). It is observed that for m > 100, the curves converge towards the 
exact time evolution, and for m=200, the DTWT curve is coincident with the exact curve. 
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FIG. 8: (Color online) Comparison of exact versus DTWT time evolution of charge (left curve) 
and spin (right curve) densities at first site and last site of biphenyl system (top) and stilbene 
molecule (bottom) whose structures are given in Fig. 4. The color coding is as follows: black 
curve with circles: exact; red curve with squares=DTWT(m=150); green curve with left trian- 
gles=DTWT(m=200); blue curve with right triangles=DTWT(m=250). 



V. REAL-TIME DYNAMICS OF SPIN-CHARGE SEPARATION IN PPP 
CHAINS 



The quantities we have studied using our time- dependent DMRG algorithm are the time 
dependence of the site charge and site spin densities. The charge and spin densities at the 
ith site of a system at time t are, respectively, 

(n,(t)) = {tPit) I {n,^ + ni_^) \ tlj{t)) (17) 

and 

(sm) = m) I (r2i.-n,_.) (18) 

where the number operators with spin a and | ip{t)) is the wave packet evolved in 

time. Although we have calculated charge (spin) density at all sites of the systems studied 
in this work, we focus on (ni(t)), {ni(t)), (s^(t)), and (s|^(t)), which suffice to investigate 
spin-charge separation in the PPP model; 1 and L correspond to the two terminal sites of 
the chain. We have done two sets of calculations for each of the systems mentioned above, 
keeping the DMEV basis sizes at 200 and 250, and we find the results converge. Here, we 
present data obtained with the smaller DMEV basis, namely, with 200 states. 

Consider the ground state of a half-filled system with particle density, 77, = (^^) = 1/2, 
A^e and A^^ are, respectively, the total number of electrons and sites in the system. The 
ground state of this neutral system being an eigenstate of the Hamiltonian, is a stationary 
state. When an electron with a definite spin is either added or removed at a site, from 
the ground state of the system, a new state forms which is no longer an eigenstate of the 
governing Hamiltonian. This is a transient (non- stationary) state, popularly known as "wave 
packet" in the literature, whose time evolution we study for a finite time. In our studies, 
we annihilate an up-spin electron from site 1 of the systems we have considered and this 
corresponds to injecting a downspin hole into the system at site 1. In the singlet {{S^^^) = 
0) ground state of the neutral system, {rii) ~ 1.0 and (s^) ~ 0.0, W i E [l,Ns], Ns = total 
number of sites in the system. The injected hole at t = is localized at the injection site 
(site 1). In the initial state (wave packet), except at the injection site, all other sites have 
charge (spin) density of 1.0 (0.0). Evolution of this non-stationary state in time, leads to 
change in the charge (spin) density distribution of the system, and temporal variation in 



these dynamical quantities is viewed as propagation of the injected hole from site 1 to site 
L. Since our systems are not connected to reservoirs, the particle number is fixed and the 
hole is reflected from the two ends of the system. Hence, time evolution profiles of charge 
(spin) density consists of a series of small and large maxima and minima. All the systems 
studied by us are homogeneous, bipartite lattices, implying that all sites are equivalent and 
the lattice possesses electron-hole symmetry at half filling. 

In the absence of electron-electron correlation, temporal variation in charge and spin 
densities of the injected hole are identical. This is because the charge and spin degrees 
of freedom of the hole propagate with the same velocity, namely, the Fermi velocity {'dp)- 
Electron-electron correlation decouples these intrinsic degrees of freedom into two separate 
elementary excitations: holon (carrying charge but no spin) and spinon (carrying spin but 
no charge). This decoupling is known as spin-charge separation and has been widely studied 
in the literature. Thus, in the presence of electron-electron correlation, the time evolution 
profiles of charge and spin densities of the system are different from each other and is 
recognized as a manifestation of spin-charge decoupling. In order to address the issue of 
spin-charge separation in the PPP model for a given topology, we have focused on two major 
extremal points in the time evolution profiles of charge and spin densities of the injected 
hole. These give us an estimate of the velocity of the charge and spin of the hole, and 
are therefore helpful in analyzing the spin-charge separation phenomena in the PPP model. 
These points correspond to the first major minima (dip) in the time evolution profiles of 
{riiit)) and (s^(t)), and {niit)) and {sj^{t)), respectively, at t 7^ 0. The time taken for the 
t ^ first minima to appear in the time evolution profiles of {ni{t)) and (sf (t)) is associated 
with the event of the charge and spin degree of freedom of the hole propagating from site 1 
to site L and returning to site 1, respectively. The time taken for this dip to appear in the 
charge (spin) density profile is designated as r^^;^ (^II)- The charge velocity deduced from 
this time is '&2l = (^) while the spin velocity is 'd^j^ = (-^). The first minima in the 
temporal profiles of {riL^t)) ((s|^(t))) appears when the charge (spin) of the injected hole 
migrates from site 1 to site L of the system. The time taken for this event is denoted as 
(r£); charge and spin velocities associated with this event are, 'd'l = (-^) and 'dl = (-^), 
respectively. 



Spin-charge separation in correlated one- dimensional systems has been studied using the 
DMRG technique by Kollath et al^ and Ulbricht and Schmitteckert^. Kollath et al. studied 
the dynamics of a wave packet obtained by introducing a particle at t = in the middle of 
a long Hubbard chain equilibrated in a spin-dependent site energy acting on the system at 
time t < and turned off at t = 0. Both for different fillings and different chain lengths, 
they observed that the charge and spin velocities were different. Ulbricht and Schmitteckert 
observed spin-charge separation in a transport simulation involving non-interacting leads and 
interacting system. The initial wave packet consists of a particle with Gaussian probability 
distribution, moving towards the interacting region, added to one of the leads. The time- 
dependent study of site spin and charge densities show a separation of spin and charge. Our 
studies are carried out on molecular systems with long-range interactions and the charge 
injection, is made at the end of the chain. 

For studying dynamics in the PPP model, we have considered regular polyene chains 

(uniform transfer integral) . This implies that all bond lengths are equal. Although polyenes 

are experimentally known to exist in dimerized form, yet our focus being on the algorithm, 

we have not considered dimerized chains in the present study. Analytical expressions for the 

velocity of the charge degree of freedom (holon), -i^/j, and spin degree of freedom (spinon), 

in the large U limit of the one- dimensional Hubbard model exist in literature^i^, 

sinf27™) 



2ti I t 

■(9/1 = 2 I t I sin(7rn); 'ds = 



2n 



(19) 



U 

where t and U are, respectively, the nearest-neighbor hopping matrix element and the on-site 
Coulomb repulsion term, and n is the particle density {n < 1). From the above expressions 
it is evinced that while -f^/i oc | t |, "i^s oc Thus, for a given value of t, increasing U is 
tantamount to decreasing the velocity of the spinon. Furthermore, as f/ — ?■ oo, — )■ 
and we reach the atomic limit. However, analysis of the charge {'dh) and spin {'dh) velocities 
of the injected hole for the one-dimensional PPP model do not exist. The PPP model is 
basically considered as a Hubbard model augmented with an additional long-range Coulomb 
repulsion term, which leads to a renormalized U in the Hubbard model. However, the physics 
of the PPP and Hubbard models are quite different; the potential energy of a configuration 
not only depends on the number of doubly occupied sites but also on the actual distribution 
of these sites. Thus, it is not possible to naively map the PPP models to effective Hubbard 
models. This makes a comparison of the charge and spin velocities of PPP model with those 




time (fs) 

FIG. 9: Time evolution profiles of (ni(i)) (left curve) and (sf(t)) (right curve) for regular PPP 
chains of length 14, 20, 30, and 40. The position of the and are indicated with arrow. 



of the Hubbard model interesting. We have carried out this comparison of the PPP results 
with those of the Hubbard model with = 2.0, 4.0, and 6.0, which we have published 
previousl)*^^. 

In Figs. 9 and 10 we show the time evolution of {ni(t)), (s^(t)) and {nL{t)), (s|(t)) for 
different chain lengths of the PPP model. Time taken for the second minimum to appear 
in the temporal profiles of {ni(t)) and (s^(t)) (Fig. 9) gives r^^. rl^^, respectively, and 
are the time taken by the charge and spin to travel from site 1 to site L and back to site 
1. Two observations can be made from the results: r^^. > '^2L chains considered. 





FIG. 10: Time evolution profiles of {nilt)) (left curve) and {sf^{t)) (right curve) for regular PPP 
chains (L = 14,20,30,40). The position of the and t£ are indicated with arrow. 



that is, the hole charge moves faster than the hole spin in the PPP model also, just as in 
the case of the Hubbard model and the velocities i)2L "^21 weakly dependent on the 
system size. Similarly, time taken for appearance of the first minimum in the time evolution 
profiles of {riiit)) and {sl{t)) (Fig. 10) provides and r£ respectively. The charge and 
spin velocities of the hole, ■^'l and ^j^, calculated from the dynamics of the Lth site as 

and ^P"^ ' respectively, agree with the velocities calculated from the r^^;^ and r|j^ values 
{^2L "^21)^ corresponding to the dynamics of the first site. We also note that finite-size 
effects are weak from the linear dependence of t^^^ on the system size and nearly system 




size independent ratio of '&/^2L (-^^S- -'--'-)• standard PPP parameters, the charge 

velocity is nearly twice the spin velocity (see Table. I). To further investigate the finite-size 
effects on correlation, we plotted the ratio of hole (spin) velocity for the 40-site chain to that 



is nearly constant, with a value close to 1.0 signifying infinite chain behavior. However, the 

spin velocity ratio shows a dip at = 4.0 and increases significantly as increases. Thus, 

it appears that the spin properties show stronger finite-size effects than charge properties. 

\ ~ liF^ ~ 1.1, and the 
J Y20J 

finite-size effects are weak. The finite-size effect in the Hiickel model arises from the kinetic- 
energy term, longer the chain, lower the kinetic energy due to delocalization, and the system 
is stable. In the Hubbard model, the finite-size effects arise due to suppression of charge 
fluctuations which is more effective in longer chains due to delocalization. However, in PPP 
models, the charge fluctuations are better accommodated due to long-range interactions and 
the kinetic term is not as dominant as non-interacting models. Therefore, we can anticipate 
weaker finite-size effects in the PPP model than either the Hiickel or Hubbard model. 

In the extreme case of correlation strength in the Hubbard model — ?■ 00, the charge 
velocity will be a constant near n ^ 0.5 and the spin velocity tends to zero [Eq. (19)]. 
The PPP model therefore appears to be in the intermediate correlation regime. In order 
to compare the PPP model with the Hubbard model, we have also carried out the time 
dependent study for different values for the Hubbard model. In Fig. 13 we compare the 
time evolution of (?t.l(^)) and {s1{t)) of the Hubbard model with = 2.0, 4.0, and 6.0, 
and the PPP model for different chain lengths. We find that for all the three correlation 
strengths of the Hubbard model considered, the charge as well as spin, move slower than 
in the PPP model. The ratio of the charge velocity to the spin velocity for the Hubbard 
model is presented in Table I. We see that for the Hubbard model with strong correlation 
strength, = 6.0, this ratio is closer to that of the PPP model. It is usually the practice 
to compute the effective in an extended range model to be |^p^, where Vi^2 is the first 
neighbor electron-electron interaction range. This simplistic interpretation of long-range 
correlations only renormalizing the U value of the Hubbard model seems to be erroneous in 
treating the djTiamics of doped holes. We note that very strong on-site correlation of > 




J versus (Fig. 12). In the case of holes it is found that this ratio 
0.0 and significantly decreases as ^ increases and beyond Ur = 4.0 it 
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FIG. 11: Variation in '7"^'^2L' ^L/2L 



represent '&2L (top) and (bottom). 



6.0 is required to reproduce the ratio of the charge and spin velocities in the PPP model. 
However, Ue// computed as for standard parameters is 1.52, underestimates severely 

the strength of correlations. It appears that the dependence of correlation on energy of the 
actual charge distribution enhances the role of correlations in the dynamics of charge and 
spin transport. 
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FIG. 12: Finite size effect of hole and spin velocities as a function of correlation strength Ratio 

-^f^ I , is plotted on the 

Y-axis. The ratios for the PPP model are quoted in the figure. 



VI. SUMMARY AND OUTLOOK 



We have developed a time- dependent DMRG scheme called DTWT technique by com- 
bining the salient features of the LXW and TST algorithms. This time-dependent DMRG 
technique is free from the drawbacks associated with both the parent methods while pos- 
sessing their strengths. Our scheme is faster and more accurate than both the LXW and 
TST methods. The TST technique targets one time step of length At and evolves the 
wave packet | ip{t)) over this time step, using four states constructed over a time interval 
At which is larger than At. The time step At is used for constructing the Hilbert space 
representing the time propagating wave packet. However, our studies have revealed that this 
usually short time-step is unsuitable for representing the Hilbert space of time-evolving sys- 
tem, the reason being. At has poor information about the future states, along the trajectory 
of time propagation of the initial wave packet. Hence, the adaptively constructed Hilbert 
space of the time evolving wave function fails to follow the evolution successfully. The LXW 
scheme targets the whole time evolution interval at every system size of the infinite DMRG 
algorithm. Hence, even though the dimension of DMEV basis (m) is fixed, the adaptively 



TABLE I: Variation in r^, r£, i?^, t^I with chain length, in the Hubbard model with ^=2.0,4.0,6.0, 
and the PPP model with t=-2A eV and [7=11.26 eV. 



Model Parameter 


L 


20 


30 


40 


^ -20 

\t\ 


" L 


6.91679985 
9.72179978 
2 024057411 
1.44006257 
1.405534352 


9.78119978 
12.8567997 
2 044738933 
1.55559707 
1.314439949 


12.3353997 
15.8069996 
3 242699951 
2.53052452 
1.281433918 


^ - 40 

1*1 


rt 
ri 


7.02239984 
12.7445997 
5 696058457 
3.138584258 
1.814849623 


10.0517998 
17.1599996 
5 969080283 
3.496503578 
1.707156921 


13.0811997 
19.0409996 
6 1 1 5647023 
4.201460096 
1.455600406 


= 60 


ri 
ri 

^ L 


7.18739984 
16.4669996 
2 782647473 
1.214550342 
2.291092742 


10.2827998 
18.9815996 
2 917493346 
1.580477970 
1.845956351 


12.9755997 
19.9319996 
3 082709156 
2.006823239 
1.536113942 


PPP 


ri 
ri 


1.81499996 
3.76859992 
11.019283989 
5.307010674 
2.076363638 


2.58059994 
5.44499988 
11.625203711 
5.509641995 
2.109974427 


3.45179992 
6.99599984 
11.588157172 
5.717553018 
2.026768643 



built Hilbert space at every system size successfully follows the time evolving wave function. 
However, the DMEV cut-off is significantly less than the number of retained target states. 
Hence the LXW procedure is capable of constructing the desired Hilbert space successfully. 
However, the time evolution of the final system is inaccurate due to the fact that number of 
target states is usually much larger than the DMRG cut-off. Furthermore, the LXW scheme 




A 



V 



0! 

'§ -0.15 
0.3 















10 15 
time (fs) 




FIG. 13: (Color online) Comparison of the time evolution profiles of (n/,(t)) (top curve) and (s|,(t)) 
(bottom curve) for regular PPP chains with those of regular Hubbard chains, for chain lengths of 20, 
30 and 40 sites. The color coding with symbols, is as follows: black curve with circles: PPP model; 
red curve with squares=Hubbard model with -^=2.0; green curve with up triangles=Hubbard 
model with -^=4.0; blue curve with down triangles=Hubbard model with -^=6.0. The position of 
T^^^ is shown with the aid of arrows: PPP model (black arrow), Hubbard model with ^=2.0 (red 
arrow), 4.0 (green arrow), and 6.0 (blue arrow). 



applied to a finite-system size, is not quasiexact. DTWT circumvents these problems ef- 
ficiently by considering a double time window (2At) for basis adaptation and single time 
window (At) which is embedded within the former time window, for evolution. The extra 
pAt steps within a double time window which are used for targeting, ensures that for every 



single time window propagation, the basis gains the information about the next single time 
window. This ensures higher accuracy of DTWT technique compared to both the LXW and 
the TST schemes. 

We have found that this method is applicable not only to chains but also to lattices with 
other topologies like rings and ring-chain systems. Using this technique we have performed 
non-equilibrium dynamics of spin-charge transport in the PPP model which harbors long- 
range Coulomb repulsion. Real-time dynamics of spin and charge transport in these systems 
is still in its infancy in the literature. In future we intend to address the effect of dimerization 
on the spin-charge separation phenomena as well as, study spin and charge transport in Y- 
j unctions. Using our DTWT procedure, we intend to study non-equilibrium transport in a 
molecule between two leads (having lead- molecule- lead geometry). 

Further improvements of the time evolution in the DTWT scheme can be carried out. 
We recognize that the time-dependent DMRG schemes have two sources of error: (i) the 
truncation error due to Hilbert-space truncation and (ii) time propagation error due to finite 
time step employed in the numerical solution of the TDSE. Besides, there is also the problem 
of stability of the numerical method for time steps larger than a critical value. We intend to 
overcome the latter two problems by using the Chebyshev polynomial-based decomposition 
of the time evolution operator. Exact propagation of the wave packet | ip{t)) using 

M 

\^{t + At)) = ^c„e-^^*(^"-^^) (20) 

n=l 

where | ipit)) = c„ | 4>n{'t)), I 4>nit)) is the nth eigenstate of i^oMRG; Ei is the ground 

state energy. En is the nth eigen-energy, M is the order of the matrix //dmrg! and c„ = 
{(pn I i^it)), has no error associated with solving the TDSE. Chebyshev expansion of the time 
evolution operator—*^ has the same advantages as exact expansion of the time evolution 
operator. The Chebyshev expansion of the state \ ip{t + At)) is given by 

I il;{t + At)) = e-™ I ijj{t)) 

M (21) 

^ ^a„T„(H) I ^(t)), 

n=0 

where T„(]H) is the nth Chebyshev polynomial of the first kind, H is the scaled Hamiltonian 



with eigenvalues ranging from [—1.0, 1.0], a„ is coefficient of T„(]H) given by 

a. = (2 - ^.0)6-^^ '""'^^^"'"'"' {-^TU^t ^-^ - ^- ). (22) 

-Emm and Emax are the minimum and maximum energies of H; is the nth order Bessel 
function of the ffist kind. The Chebyshev expansion of the evolution operator can be evalu- 
ated up to machine accuracy and there will be no error associated with this time propagation 
scheme for any arbitrary step-size. The use of Chebyshev polynomial expansion scheme is 
expected to render the DTWT technique free from any time step error besides decreasing 
the computational time. 
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